Evaporation of Lennard-Jones Fluids 
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Evaporation and condensation at a liquid/vapor interface are ubiquitous interphase mass and 
energy transfer phenomena that are still not well understood. We have carried out large scale 
molecular dynamics simulations of Lennard-Jones (LJ) fluids composed of monomers, dimers, or 
trimers to investigate these processes with molecular detail. For LJ monomers in contact with a 
vacuum, the evaporation rate is found to be very high with significant evaporative cooling and an 
accompanying density gradient in the liquid domain near the liquid/vapor interface. Increasing the 
chain length to just dimers significantly reduces the evaporation rate. We confirm that mechanical 
equilibrium plays a key role in determining the evaporation rate and the density and temperature 
profiles across the liquid/vapor interface. The velocity distributions of evaporated molecules and 
the evaporation and condensation coefficients are measured and compared to the predictions of an 
existing model based on kinetic theory of gases. Our results indicate that for both monatomic and 
polyatomic molecules, the evaporation and condensation coefficients are equal when systems are not 
far from equilibrium and smaller than one, and decrease with increasing temperature. For the same 
reduced temperature T/T c , where T c is the critical temperature, these two coefficients are higher 
for LJ dimers and trimers than for monomers, in contrast to the traditional viewpoint that they are 
close to unity for monatomic molecules and decrease for polyatomic molecules. Furthermore, data 
for the two coefficients collapse onto a master curve when plotted against a translational length 
ratio between the liquid and vapor phase. 



I. INTRODUCTION 

The inverse processes of evaporation and condensation 
are of fundamental importance in natural phenomena 
and engineering applications. In both processes, heat 
and mass transfer between liquid and vapor phases^ The 
key physical quantities to determine are the interphase 
mass and energy transfer rates. There have been a num- 
ber of theoretical analyses, using either kinetic theory^— 
or nonequilibrium thermodynamics,— ~— to predict these 
rates as well as the density, temperature, and pressure 
profiles in the vapor phase. In the framework of ki- 
netic theory of gases, the problem was typically formu- 
lated based on the Boltzmann-Bhatnagar-Gross-Krook- 
Welander (BBGKW) equation for the vapor, with the 
liquid surface temperature and the temperature, pres- 
sure, and velocity of vapor far away from the inter- 
face as boundary conditions. To obtain solutions of the 
BBGKW equation, one further assumption of the liq- 
uid/vapor interface was usually made, i.e., all molecules 
approaching the interface completely condense into the 
liquid phase, while molecules evaporated from the liq- 
uid surface have a Maxwell-Boltzmann (MB) distribu- 
tion corresponding to the saturated vapor at the liquid 
temperature. This is equivalent to assuming that the 
evaporation and condensation coefficients are unity since 
the evaporation (condensation) coefficient is defined as 
the ratio of an experimental evaporation (condensation) 
rate to a theoretical maximum rate given by the Hertz- 
Knudscn (HK) equation, which is exactly the rate cor- 
responding to a MB distribution. Though the evapora- 
tion and condensation coefficients are closely related and 
have the same value for an interface in equilibrium, they 
could be substantially smaller than unity and have differ- 
ent values for nonequilibrium interfaces. The above as- 



sumption of the boundary condition at the liquid/ vapor 
interface might not be realistic and recently some theo- 
retical work emerged attempting to replace it with more 
physical ones^~— 

The kinetic theory of evaporation and condensation 
was challenged in a recent experiment of Fang and 
WardJ^ They measured the temperature profile to within 
one mean free path (~ 19/im) of the interface of an 
evaporating liquid and found a discontinuity in temper- 
ature across the interface that was much larger in mag- 
nitude and in the opposite direction to that predicted by 
the kinetic theory or nonequilibrium thermodynamics. 
This disagreement led Fang and Ward to conclude that 
the boundary conditions traditionally assumed for the 
BBGKW equation were unphysical and they developed 
a statistical rate theory of evaporation flux to explain 
their experimental observations. 20 ' 21 

To elucidate the physics of liquid/ vapor interface dur- 
ing evaporation and condensation, more detailed mea- 
surements are still needed at even smaller scales and with 
more refined resolutions. However, the liquid/vapor in- 
terface is difficult to probe experimentally and available 
data can only provide limited information of the micro- 
scopic detail of an evaporating interface. In the past sev- 
eral decades, although many measurements have been re- 
ported on the evaporation and condensation coefficients 
of various materials, it still remains unclear how these re- 
sults can be related to the molecular processes occurring 
at the liquid/vapor interface. Furthermore, reported val- 
ues are often scattered in such a large range even for the 
same material which further complicates their theoretical 
interpretation. For example, the reported measurements 
of the evaporation coefficient of water range from 0.01 to 
1.0 ) 22 i 23 which exemplifies the difficulty to obtain accu- 
rate values of these coefficients, let alone the molecular 
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mechanism of evaporation and condensation. 

The development of computer simulation techniques, 
particularly the molecular dynamics (MD) method, has 
enabled a number of studies of the evaporation and 
condensation processes at the molecular scale ) 16 i 17 ' 24 ~ — 
which have advanced substantially our understanding 
of these interphase mass and energy transfer phenom- 
ena. Matsumoto et al. studied the liquid-vapor inter- 
faces of argon, water, and methanol, and obtained the 
evaporation and condensation coefficients that agreed 
reasonably well with experimental values^ir— Tsuruta 
et al. studied the condensation coefficient of Lennard- 
Jones (LJ) fluids by injecting test molecules to bombard 
the liquid/vapor interface and measured the reflecting 
probability, from which the condensation coefficient was 
deducedi22r— Their results confirmed that the evapora- 
tion and condensation coefficients are equal in equilib- 
rium systems, but revealed that the condensation prob- 
ability generally depends on the normal component of 
the kinetic energy of incident molecules. This is in con- 
trast to the common assumption that the condensation 
probability is constant for vapor molecules that hit the 
liquid surface. They also measured the velocity distribu- 
tions of the evaporated and reflected molecules at the liq- 
uid/vapor interface. Inspired by the MD results, Tsuruta 
et al. developed an expression for the velocity-dependent 
condensation probability and later justified this expres- 
sion using the transition state theor y 30 ' 31 which was used 
earlier to estimate the condensation coefficient by Fu- 
jikawa and Maerefati^ Anisimov et al. also used MD to 
investigate the evaporation of LJ fluids and the proper- 
ties of liquid/vapor interface, and extended their studies 
to the case of high-rate evaporation ^ 3 ' 34 

Recently, Rosjorde et al. used MD simulations to study 
L J / spline fluids and found evidence for the hypothesis of 
local equilibrium at a liquid/vapor interface i 35 ' 36 They 
also measured transfer coefficients of the mass and en- 
ergy fluxes and found that they agreed with kinetic the- 
ory from the triple point to about halfway to the critical 
point. They suspected that the disagreement between ki- 
netic theory predictions and experimental results of Fang 
and Ward was due to the fact that the theory dealt with 
monatomic fluids, while in the experiment polyatomic 
fluids were used. Meland et al. also studied LJ/spline 
fluids and compared their MD results with gas-kinetic 
calculations^ They found that the evaporation and con- 
densation coefficients are not equal outside equilibrium 
and there is a significant drift velocity in the distribu- 
tion function at the interphase for both net evaporation 
and net condensation. Ishiyama et al. used MD simula- 
tions of LJ fluids to check the validity of kinetic bound- 
ary condition for the BBGKW equation and found that 
the condensation coefficient is close to unity below the 
triple-point temperature and decreases gradually as the 
temperature rises.— 

More recently, Holyst and Litniewski studied the evap- 
oration of nanodroplets and demonstrated that the evap- 
oration process is limited by the heat transfer and energy 



balance condition^ This finding challenges the basic as- 
sumption of kinetic theory that the evaporation flux is 
determined by the diffusion of mass in the vapor phase. 
In another study, they simulated a LJ liquid film evapo- 
rating into a vacuum and measured the density, tempera- 
ture, and pressure profiles^ Holyst and Litniewski found 
that mechanical equilibrium is established very quickly 
and derived an expression for the mass flux that described 
their simulation results much better than the frequently 
used HK formula. 

In these previous studies, only simple LJ fluids com- 
posed of monomers were used. However, it is well known 
that the vapor pressure of a LJ fluid is much higher than 
those of real liquids. Though a high vapor pressure im- 
plies that LJ liquids are relatively easy to evaporate, it 
does not necessarily mean that the evaporation coeffi- 
cient is close to unity because the corresponding theoret- 
ical maximum flux is also large. In fact, some previous 
simulations found that at a moderate temperature the 
evaporation and condensation coefficients for simple LJ 
liquids are around 0.8^^^. 

In this paper, we first simulate the evaporation of LJ 
monomers and then extend MD simulations to the evapo- 
ration of LJ fluids composed of dimers and trimers. This 
allows us to examine the effect of molecular composi- 
tion on the evaporation rate and the evaporation and 
condensation coefficients. The paper is organized as fol- 
lows. In Sec. II, the simulation methodology and the 
procedure for creating the liquid/ vapor interface and re- 
moving atoms to implement controlled evaporation are 
introduced. Results on phase diagrams and surface ten- 
sions of dimers and trimers are presented in Sec. III. In 
Sec. IV results on density and temperature profiles are 
presented. The measured evaporation rates are found 
to agree well with the modified HK expression derived 
by Holyst and Litniewski 4i The role of mechanical equi- 
librium during evaporation is also examined. Then in 
Sec. V we measure the velocity distributions of the evap- 
orated and condensed molecules and compare them to 
predictions of a kinetic model of evaporation based on 
the transition state theory of Tsuruta et al. The evapo- 
ration and condensation coefficients are also determined. 
A brief summary and conclusions are included in Sec. VI. 

II. SIMULATION METHODOLOGY 

We carried out large scale MD simulations of three 
simple liquids. In all three cases the interaction between 
atoms is described by the standard LJ 12-6 interaction 

U(r) = 4e [(a/r) 12 - {a/rf - (a/r c ) 12 + (<r/r c ) 6 } , 

(1) 

where r is the distance between two atoms, e is the unit 
of energy, and a is the diameter of the atom. The interac- 
tion is truncated at r c = 2.5a. While commonly used to 
simulate liquids, the LJ model has one serious drawback 
when modeling the evaporation process. Namely it has 
a very high vapor pressure. This results in a very large 
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vapor density which does not correspond well to most liq- 
uids. One simple way to make the model more realistic 
but still retain the simplicity of the LJ interaction is to 
model the solvent as short chains of LJ atoms. In this pa- 
per we show that by modeling dimers that consist of two 
atoms and trimers consisting of a linear chain of three 
atoms, we obtain a more reasonable low vapor density. 
In dimer and trimer molecules, the bonded atoms are 
connected by an additional finitely extensible nonlinear 
elastic (FENE) potential, 



f/ FE N E (r) = -ifc-Rgln [1 - (r/R ) 2 



(2) 



where k — "SOk^T/a 2 , T is the temperature, /cb is the 
Boltzmann constant, and i?o = 1.5(7 

All MD simulations were performed using the 
LAMMPS simulation packaged Each system started 
with TV = 3,000,000 atoms in a parallelepiped of di- 



mensions L x x L y x L z , where L a 



L y = 200fi. The 



simulation cell was periodic in the x and y directions and 
had upper and lower confining walls in the z direction as 
illustrated in Fig. [T] Both upper and lower z- walls inter- 
acted with the monomers with a LJ 9-3 potential, which 
depends only on the distance z from the wall, 



U(z) = e v 



^(a/zf {a/zf - ^/z c f + (a/z c f 



(3) 

where e w = 2e. For the lower wall the interaction was 
truncated at z c = 2.5a, while at the upper wall was 
purely repulsive with z c — 0.72cr. Each liquid film was 
constructed by placing N atoms randomly in the simula- 
tion cell with L z chosen so that the system was near its 
bulk liquid density. Overlaps were removed by running 
a short simulation with the "nve/limit" fix in LAMMPS 
turned on42 This limits the maximum displacement of an 
atom per time step and is a very efficient way to remove 
overlaps between atoms. The system was then equili- 
brated at pressure P = by adjusting the position of 
the upper wall. After the equilibration, the position of 
the upper wall was increased by at least 50er to form a 
liquid/ vapor interface as shown in Fig. Q] The system 
was then allowed to come to equilibrium with its vapor 
before the evaporation process was initiated. 

The equations of motion were integrated using a 
velocity- Verlet algorithm with a time step St = 0.005r 
for the dimers and trimers and O.Olr for the monomers, 
where r = aim/t) 1 ^ 2 and m is the monomer mass. How- 
ever, in simulations to investigate the mechanical equi- 
librium during evaporation, the time step was reduced 
to 0.001t in order to get the small pressures accurate in 
the liquid phase. During the equilibration, T was held 
constant by a Langevin thermostat weakly coupled to 
all atoms with a damping constant T — O.lr -1 . Once 
the liquid/vapor interface was equilibrated, the Langevin 
thermostat was removed except for those atoms within 
15a of the lower boundary at z = 0. We refer to the liq- 
uid temperature in this region as the bulk temperature 




thermostated 



FIG. 1: (Color online) Liquid/vapor equilibrium at T = 
0.9e/fcB for LJ monomers (left) and dimers (right). The vapor 
density is clearly lower for dimers. 



Tft. The system was then equilibrated for an additional 
20, 000 to 100, 000 time steps before the evaporation pro- 
cess was initiated. 

To model evaporation, a deletion zone of thickness 20er 
approximately 50cr above the starting liquid/ vapor inter- 
face was defined where atoms were removed at a specified 
rate. For dimer and trimer systems, once one atom from a 
molecule was removed, the entire molecule was removed. 
In this paper, the evaporation rate is defined as the num- 
ber of atoms removed per unit time and area. There- 
fore, one removed dimer (trimer) molecule contributes 
two (three) evaporated atoms. In simulations with a con- 
trolled evaporation rate, n atoms in the deletion zone 
were removed every N t time steps. Typically, n = 10 
and Nt — 100 or 1000 for monomers and dimers. Since 
each trimer molecule contains 3 atoms, setting n = 9 
and N t = 90 or 900 for trimer systems ensured the same 
evaporation rates as for the other two systems. By re- 
moving all atoms that entered the deletion zone, a system 
in contact with a vacuum was effectively modeled. Since 
only those atoms within 15a of the lower wall were cou- 
pled to the thermostat during the evaporation process, 
the thermostat did not affect the evaporation process at 
the liquid/ vapor interface. To avoid any finite size effects 
we only analyzed evaporation data for films thicknesses 
larger than ~ 50a. 

In our simulations, the density and temperature pro- 
files in the simulation cell were measured. For liq- 
uid/vapor equilibrium cases, time-averages of the den- 
sity and surface tension were carried out to determine 
the liquid/vapor phase diagrams. For systems undergo- 
ing evaporation, instantaneous density and temperature 
distributions were calculated. Since these systems are 
translationally invariant in the x-y plane, both quanti- 
ties only depend on z and results presented here were 
averaged over x and y directions. 

The meaning and definition of temperature in non- 
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equilibrium systems still have many ambiguities 42 In 
this paper, we measured a local temperature T(z) in two 
ways. In one, T{z) is taken as the mean kinetic energy 
of atoms locating in the spacial region from z — Az to 
z + Az, where Az is typically 0.5a. The corresponding 
expression is 



T(z) 



z+A 

3JVjfcn ^ 



(4) 



Az 



where Nt is the number of atoms in the region and v 
is the atomic velocity. Another definition of the local 
temperature is based on velocities relative to the possi- 
ble advective motion induced by evaporation and can be 
written as 



T r (z) 



z+Az 



z-A: 



(v - vf 



(5) 



where v is the mean atomic velocity in the spatial region 
from z — Az to z + Az. Clearly, the two definitions give 
identical results in equilibrium since v = there. How- 
ever the results are generally different in nonequilibrium 
states. This will be discussed further in Sec. IV. 

III. LIQUID/VAPOR EQUILIBRIUM 

For the system composed of LJ monomers, the liq- 
uid/vapor phase diagram is well known4^~— However, 
this is not the case for dimers and trimers. Therefore 
for these two systems we also simulated the liquid/vapor 
coexistence by removing the walls in the z direction and 
replacing them with periodic boundary conditions, which 
resulted in two liquid/ vapor interfaces in the simulation 
cell. The liquid and vapor densities ph and pv respec- 
tively along the coexistence curve for all three systems 
are shown in Fig. [2] The liquid/vapor critical temper- 
ature T c for each system is determined from fitting the 
measured densities to 



Pl + Pv = a — bT, 
p L -pv = A(l-T/T c f, 



(6) 



where a, 6, and A are fitting parameters 42. The criti- 
cal exponent j3' — 0.318 is for Ising systems that are 
in the same universality class as simple fluids; thus it 
is fixed to this value in the fitting. The best fits using 
Eq. (J6j) give the critical temperature T c = 1.085, 1.475, 
and 1.720e/fcB, and the critical density p c — 0.316, 0.299, 
and 0.294m/cr 3 for the monomer, dimer, and trimer sys- 
tems, respectively. 

The surface tension of the liquid/ vapor interface was 
determined by measuring the stress tensor. Since the 
two interfaces are parallel to the x-y plane, the surface 
tension 7 is given by the Kirkwood-Buff formula^ 



7 



3 it) Z \P"i Z ) ~ {Pxx(z) +Pyy{z))/2]dz (7) 
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FIG. 2: (Color online ) Liquid/vapor phase diagram for LJ 
monomers (triangles), dimers (circles), and trimers (squares). 
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FIG. 3: (Color online) Surface tension 7 as a function of tem- 
perature for LJ monomers (triangles), dimers (circles) and 
trimers (squares). The results for monomers are taken from 
Ref . [H 



where p xx {z), p yy {z), and p zz (z) are the three diagonal 
components of the stress tensor and P x , P y , and P z are 
the spatially averaged pressure in each direction. The 
factor 1/2 before the integral comes from the fact that 
there are two interfaces. Far from the liquid/vapor inter- 
face, both the liquid and vapor phases are homogeneous 
and isotropic, and all diagonal components of the stress 
tensor are the same and equal to the hydrostatic pres- 
sure. Therefore in these regions the integrand in Eq. Q 
is zero and does not contribute to the integral. Near the 
interface whose normal is along the z-direction, p xx (z) 
and Pyy(z) are less than p zz {z), leading to a difference in 
P x , P y , and P z . The net outcome is a nonzero surface 
tension of the liquid/vapor interface. Results for 7 are 
shown in Fig. [3] for LJ monomers, dimers, and trimers. 
As expected, 7 drops as temperature is raised and ap- 
proaches zero as T c is reached. When the temperature T 
is well below T c , the surface tension 7 roughly decreases 
linearly with T . However, when T c is approached the 
rate of reduction of 7 becomes smaller. It is expected 
7 ~ (T c — T) u with v as a critical exponent. Results in 
Fig. [3] are consistent with v > lM- 
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Figure [3] also shows that the surface tension is smaller 
for monomers and larger for dirners and trimers. This 
implies that the cohesion between molecules becomes 
stronger as the chain length gets longer. Accordingly, at 
the same temperature the vapor densities of dimers and 
trimers are lower than that of monomers, as illustrated 
in Fig.Q] Note that dimers and trimers considered in this 
paper are very flexible objects. It is not surprising that 
longer chains lead to stronger attractive interactions be- 
tween molecules because there are more interacting sites 
and the counteractive effect of rotational degrees of free- 
dom is suppressed. We will see later that the stronger 
cohesion in dimers and trimers also affects their evapo- 
ration and condensation coefficients. 
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OF LENNARD-JONES 



IV. EVAPORATION 
FLUIDS 

A. Monomers 



When placed in contact with a vacuum, the evapora- 
tion of the LJ monomer system occurs very rapidly, as 
shown in Fig. |4ja) . The vapor is quickly depleted and the 
vapor density drops by an order of magnitude from its 
equilibrium value. On the other hand, the liquid density 
increases by about 20% from its equilibrium value near 
the interface^ as there is significant evaporative cooling 
of the liquid film [Figs. SJb) and0Jc)], no matter the 
temperature is measured either as T(z) in Eq. (@| or as 
T r (z) in Eq. <j5j) . Note that there is essentially no ad- 
vective motion in the liquid region, and thus v = and 
T{z) — T r {z) there. However, since the vapor is very di- 
lute when in contact with a vacuum, all molecules leaving 
the liquid phase almost move freely toward the deletion 
zone and get removed. As a consequence of the absence 
of collision, the vapor region is far fro m local e quilibrium 
and the mean velocity v is of order ^Jk^T '/m. Therefore 
T{z) is very different from T r {z) in the vapor phase, as 
shown in Figs.|4jb) and|4jc). 

Though T(z) and T r (z) are quite different in the vapor 
phase when in contact with a vacuum, they are the same 
in the liquid region in all cases. They are also very close in 
both the liquid and vapor regions in our other simulations 
described later where the evaporation rate is controlled 
and the vapor density during evaporation remains com- 
parable to its equilibrium value. From an experimental 
perspective, it is the average kinetic energy that is mea- 
sured when a thermocouple is placed in the vapor phase 
during evaporation. With these considerations, we refer 
to T(z) as temperature in this paper and report its mea- 
surements hereafter. Holyst and Litniewskiil called this 
temperature the pseudotemperature and concluded that 
it is very useful in describing the evaporation kinetics. 

Figure S] shows that both the density and temperature 
gradients increase as the evaporation process continues. 
After a short time, the temperature has decreased across 
the interface from its initial value of 0.9e/fcB to around 
0.7e/fcB, which is very close to the triple point of this 
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FIG. 4: (Color Online) (a) Density, (b) temperature [defined 
as the average kinetic energy, Eq. (j4}], and (c) temperature 
[with mean velocity subtracted, Eq. ([5])] for a LJ monomer 
fluid at Tb = 0.9e/fcB in contact with a vacuum. From right to 
left the profiles are plotted every 2000r since the evaporation 
process was started at t = (the rightmost curve). 



system.— Similar results were found for T& = 0.8e/fcB 
and l.Oe/kB, where the temperature at the liquid/vapor 
interface also drops to ~ 0.7e/fcB- The large increase 
in density and decrease in temperature near the inter- 
face reported here for LJ monomers are much stronger 
than those of most common fluids such as water. In this 
sense, the particular feature of LJ model makes the effect 
of evaporation more dramatic. The density enhancement 
near an evaporating interface was also seen in the simu- 
lations of Holyst and Litniewski.— 

An evaporation rate je is defined as the number of 
atoms removed in the deletion zone per unit time and 
area. It is related to the total number Ne of removed 
atoms through 



JE 



I dN E 

LrpL,, dt 



(8) 



Results for j E and N E are shown in Fig. [5] for LJ 
monomers at three temperatures. As expected, je is 
larger at higher since evaporation is a thermally acti- 
vated process and occurs more rapidly at higher temper- 
ature. The rate je also shows strong time dependence. 
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FIG. 5: (Color Online) Evaporation rate ]e vs. time for 
a monomer system in contact with a vacuum at Ti = 0.8 
(solid line), 0.9 (dashed line), and 1.0e/fcB (dotted line). The 
predicted evaporation flux j m /m based on Eq. ^ is shown as 
circles. The inset shows the total number of removed atoms 
Ne vs. time. 

In all cases, je initially has a high value since the LJ fluid 
has a high vapor density and there are plenty of vapor 
molecules in the deletion zone. As the vapor is rapidly 
depleted, je drops significantly over time and eventually 
reaches a plateau. The final reduction factor, defined as 
the ratio between the value of Je at t = when the evapo- 
ration was initiated and that at very large t, is about 2 for 
Tb = 0.8e/A;B and increases to about 4 for = 0.9e/fcB 
and 1.0e/k B . 

Holyst and Litniewski demonstrated that during evap- 
oration of a liquid film into a vacuum the momentum 
flux, j p = n va _pM (u 2 ) in the vapor phase far from the liq- 
uid/vapor interface is equal to the pressure in the liquid 
film, piiq^i Here n vap is the vapor number density, M 
is the molecular mass, and (u 2 ) is the mean squared z 
component of the molecular velocity. From this observa- 
tion they proposed an equation for the mass flux during 
evaporation: 

(U z ) 

jrn = n vap M(u z ) = piiq-r-^-, (9) 

\ u z) 

where (u z ) is the mean z component of the molecular 
velocity. Note that j m has the unit of mjE- In our simu- 
lations, pii q , (w 2 ), and {u 2 z ) are averaged in a thin region 
from Zint(i) + z\ to Zi nt (t) + z 2 , where z- lnt (t) denotes the 
location of interface at time t. Typically, z\ = — 30<x and 
z 2 = —lOcr for piiqj and Z\ = 20a and z 2 = 40cr for (u z ) 
and (u 2 ). However, results are not sensitive to the ex- 
act location of these regions and varies by a few percent 
at most as long as they are far from the liquid/ vapor 
interface. The calculated mass flux j m from the above 
equation at Tf, = 0.9e//cB is shown in Fig. [3] Clearly, the 
measured je agrees with j m /m at most times, except at 
the beginning stage of evaporation when j m /m tends to 
be a few percent smaller than je ■ 

We also ran simulations in which je was controlled by 
limiting the number of atoms removed from the deletion 



zone during a given simulation period. Results for the 
density and temperature profiles are shown in Fig. [5] for 
two rates that are smaller than the final evaporation rate 
in the case of contacting with a vacuum, the plateau in 
Fig. by a factor around 4 and 16, respectively. In 
Figs. HJa) and[SJb), the density enhancement and tem- 
perature drop are still clearly visible, but the magnitudes 
are greatly reduced compared to those in Fig. SJ The 
vapor density in this case is also reduced from the equi- 
librium value but remains finite. For a very slow evapo- 
ration rate [Figs.^c) and|5Jd)], these two effects almost 
disappear and the vapor density is close to the equilib- 
rium density. 
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FIG. 6: (Color Online) Density [(a) and (c)] and tempera- 
ture [(b) and (d)] for a LJ monomer fluid at Tt, = 0.9e/fcB 
evaporating at a rate je = 1.0 x 10~ 3 t~ 1 <j~ 2 [(a) and (b)] or 
j E = 2.5 x 10~ 4 T _1 cr~ 2 [(c) and (d)]. From right to left the 
profiles are plotted every 8000t [(a) and (b)] or 24, OOOr [(c) 
and (d)] since the evaporation process was started at t = 
(the rightmost curve). The inset in (c) shows a blow up of 
the density in the interfacial region. 

In addition to the evaporative cooling, the temperature 
profiles shown here [Figs. HJb), |^b) and Efd)] indicate 
that the temperature of vapor increases with the distance 
from the liquid/ vapor interface. This can be understood 
by a simple argument, molecules with higher kinetic en- 
ergies evaporate faster and more energetic molecules ac- 
cumulate near the deletion zone, leading to a higher tem- 
perature at this end and a density gradient in the vapor 
phase. Therefore, if the vapor temperature is measured 
at a certain distance away from the liquid/vapor inter- 
face, it will be higher than the actual temperature at 
the interface. This phenomenon is qualitatively consis- 
tent with the experimental measurement by Fang and 
Ward. 19 They measured the temperature as close as one 
mean free path of the interface of an evaporating liquid 
and the results indicated that the vapor temperature is 
indeed greater than that in the liquid phase at the inter- 
face. They also concluded this is due to more energetic 
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molecules that are likely to evaporate first. Note that 
the experimental condition in Ref. [l9| is close to the sit- 
uation where the evaporation rate is controlled here. In 
this case, the vapor density is generally large enough to 
make local thermal equilibrium approximately true. Our 
simulations confirm ed that the mean velocity v is much 
less than y/k-aT/m, which leads to T(z) ~ T r (z). Since 
the quantity recorded by the thermocouple is related to 
the average kinetic energy of vapor molecules, the qual- 
itative comparison made here between simulations and 
experiments should be reasonable. 

B. Dimers and Trimers 

While the LJ potential still serves as a reasonable ap- 
proximation of inter-molecular interactions, most liquids 
of interest are composed of molecules not single atoms. 
One simple extension is to consider molecules of two or 
more LJ monomers bound together. In the present simu- 
lations, the bonded interaction is realized through the in- 
troduction of a FENE potential [Eq. ([IJ] between bonded 
monomers. The nonbonded interatomic interactions are 
still given by the LJ potential [Eq. (JTJ]. 
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FIG. 7: (Color Online) Density [(a) and (c)] and temperature 
[(b) and (d)] for the LJ dimer fluid at Ti, = 1.2e/fcB in contact 
with a vacuum [(a) and (b)] or evaporating at a controlled rate 
j E = 5.0 x lCTV -1 ff -a [(c) and (d)]. From right to left the 
profiles are plotted every 5000Y [(a) and (b)] and 100, OOOr [(c) 
and (d)] since the evaporation process was started at t = 
(the rightmost curve). The inset in (c) blows up the density 
profile in the liquid and interfacial region. Temperature in 
the vapor region is not included because the density of vapor 
is small and the data are too noisy to indicate a clear trend. 

Results on the density and temperature profiles for the 
LJ dimer system in contact with a vacuum or evaporating 
at a controlled evaporation rate are shown in Fig. [7] for 
Tf, = 1.2e/fce. Those for the trimer system at the same 
bulk temperature but at two controlled rates are shown 
in Fig. [9J Results on the evaporation rate je at various 



temperatures when in contact with a vacuum are shown 
in Fig. [5] (dimer) and Fig. [TO] (trimer). 

The phenomenology of evaporation into a vacuum for 
dimers [Fig. [7fa) and EJb)] and trimers (not shown) 
is similar to that of monomers, including the density 
enhancement and the evaporative cooling near the liq- 
uid/vapor interface. However, the relative magnitudes 
of density increases and temperature drops are smaller 
for dimers and even smaller for trimers. This is consis- 
tent with the observation that cohesion gets stronger and 
evaporation slows down in dimer and trimer fluids. 
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FIG. 8: (Color Online) Evaporation rate js vs. time for 
the LJ dimer system in contact with a vacuum at Tt = 0.9, 
1.0, 1.1, and 1.2e/fcB (from bottom to top). The predicted 
evaporation flux jm/m based on Eq. @ is shown as circles. 
The inset shows the total number of removed atoms Ne vs. 
time. 

More quantitative information of evaporation into a 
vacuum is obtained by measuring the evaporation rate 
jE- It generally decreases with time and only approaches 
a constant value after a certain interval. Compared with 
monomers the reduction factor is smaller for dimers. For 
the trimer system the time-dependence is weak except 
near the critical temperature T c . This is understandable 
because the vapor density decreases as the chain length 
increases. As expected, je is higher at higher T^. At the 
same bulk temperature Xb, je is largest for the monomer 
system and smallest for the trimer. For example, at T& = 
0.9e/A;B it is around 2 x 10 _3 t _ V -2 for monomers, 4 x 
10~ 4 T <7 -2 for dimers, and 4 x 10 _5 T~ 1 er~ 2 for trimers. 
Thus increasing the number of atoms in a molecule only 
from 1 to 3 results in a decrease of more than a factor of 
50 in ]e ■ 

The results for j m / m from Eq. ([9]) for dimers at TJ, = 
1.2e/fce and for trimers at Ty, = 1.2e/k^ and 1.5e/fcB are 
shown in Figs. [5] and [TUJ respectively. In all cases, the 
agreement between j m /m and je is quite satisfactory, 
indicating that Eq. (J9j) is valid not only for monatomic 
liquids, but also for molecular liquids such as dimers and 
trimers. We also note that in evaluating (w 2 ) in Eq. ©, 
molecular velocities should be used. 
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FIG. 9: (Color Online) Density [(a) and (c)] and tempera- 
ture [(b) and (d)] for the LJ trimer fluid at Tj, = 1.2e/feB 
evaporating at a rate ]e = 5.0 x 10 _4 t _1 <7 -2 [(a) and (b)] or 
je = 5.0 x W 5 T^ 1 a' 2 [(c) and (d)]. From right to left the 
profiles are plotted every 60, OOOr in both cases since the evap- 
oration process was started at t = (the rightmost curve). 
The insets in (a) and (c) show a blow up of density in the 
liquid and interfacial region. 
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FIG. 10: (Color Online) Evaporation rate je vs. time for the 
LJ trimer system in contact with a vacuum at Tb = 0.9, 1.2, 
and 1.5e/fea (lines from bottom to top). Open squares and 
circles show the predicted evaporation flux j m /rn from Eq. Q 
at Tb — 1.2 and 1.5e/feB, respectively. The inset shows the 
total number of removed atoms Ne vs. time. 



The magnitude of the density enhancement and tem- 
perature drop near the interface also depends on the 
evaporation rate. As shown in Fig. [7f c) and (d), for the 
dimer system at Tb = 1.2e/fcB they are not observed when 
]e is reduced to 5.0 x 10 _5 T _1 er~ 2 , which is much smaller 
than the plateau value of je of this system in contact 
with a vacuum. For the trimer system under the same 
temperature, the density and temperature gradients are 
visible at j E = 5.0 x lO^T^a" 2 [Fig.^a) and^b)], but 
disappear when je is reduced further by a factor of 10 to 



j E = 5.0 x lO-'V-V- 2 [Fig. Etc) andEtd)]. Note that 
the plateau value of j E for trimers at this temperature 
when evaporating into a vacuum is 7.7 x 10 _4 t _1 (7~ 2 . 
Thus it is reasonable to see density and temperature gra- 
dients at js — 5.0 x 10 _4 r _1 cr -2 , but not at much slower 
evaporation rates. 

C. Stagnation Pressure During Evaporation 

During evaporation, the vapor escaping into a vac- 
uum is far from equilibrium and a local temperature and 
pressure defined with respect to a thermodynamic equi- 
librium cannot be used for its description. Holyst and 
Litnicwski showed for LJ monomers that global quanti- 
ties, such as total kinetic energy, total mass flux, and 
total momentum flux, can be well defined and provide 
useful information about the physical state of the vapor 
phased Particularly, a stagnation pressure tensor can be 
defined as 

Pa ? = T7 \J2 "hm^iP - E £ S^- ' ( 10 ) 

V V 3 i>0 3 ° Tij0 J 

where i and j index all atoms inside the volume V, rrij 
and Vj are the atomic mass and velocity, and </>y are 
the interatomic distance and potential, respectively, and 
Greek indices indicate components along the x, y, and z 
directions. Note that this pressure tensor is not defined 
relative to a local equilibrium. Holyst and Litniewski 
demonstrated that for the stagnation pressure the p zz 
component in the vapor phase is equal to the liquid pres- 
sure, piiq, far from the interface.— 

We measured the stagnation pressure for Lennard- 
Jones monomers, dimers, and trimers. Results are shown 
in Fig. [TT] Indeed, for all systems at all times the p zz 
component in the vapor phase (solid lines in the main 
panel of Fig. [TTj) is equal to the pressure in the liquid 
film (symbols) . Data in Fig. [TT] are for liquid films evap- 
orating into a vacuum. We confirmed that the equality 
is also true for systems evaporating at controlled rates. 

The inset of Fig. [TTJ shows the stagnation pressure 
along the normal direction of the liquid/ vapor interface. 
In the liquid phase, p zz = p xx = p yy = py lq . However, 
in the vapor phase, p zz ^ p xx — p yy , and the difference 
grows as the distance away from the interface. The inset 
also shows that the momentum flux, j p = n vap M(u 2 ), 
approaches p zz in the vapor phase far from the interface. 
This is not surprising for monatomic systems since in 
this case j p is exactly the first term in the expression for 
p zz in Eq. (jlpp . The second term in p zz represents the 
contribution from interatomic interactions. For systems 
evaporating into a vacuum, the vapor becomes increas- 
ingly dilute away from the liquid/ vapor interface so that 
the contribution to p zz from interactions between two 
monomers essentially vanishes and the momentum flux 
dominates. For molecular systems, j p is not identical to 
the first term in p zz as expressed in Eq. (|10p . which is 
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FIG. 11: (Color Online) Stagnation pressure p zz averaged 
in the liquid phase (symbols) and vapor phase (lines) for: 
monomers at Tb = 0.9e/fcB (the middle line and 0)i dimers 
at Tt — 1.2e/fcs (the bottom line and □); and trimers at 
Tb = 1.5e/fcB (the top line and v)- The trimer data are 
shifted upward by 0.002e/<r 3 and right by 200r for clarity. 
The inset shows the stagnation pressure p as a function of 
distance along the normal direction of the liquid/vapor inter- 
face for trimers at Tb = 1.5e/fcB- The pressure in the liquid 
phase, pn q , is shown as squares (□), and that in the vapor 
phase, Pvap, is shown as the dashed line (the p zz component) 
and the dotted line [the (p xx +p yy )/2 component]. The solid 
line represents the momentum flux j p = n va , p M{u 2 z ). The in- 
terfacial region is denoted as a shadowed strip. Data are for 
systems evaporating into a vacuum. 



based on atomic velocities not molecular velocities aver- 
aged from those of its constitute atoms. However, p zz 
can be also defined in terms of molecular velocities and 
the interactions between molecules instead of atoms^ 
Then molecular density becomes dilute during evapora- 
tion into a vacuum and the same argument above applies 
again. The result is that j p still dominates over the in- 
teraction contribution in p zz in molecules escaping into 
a vacuum. This relation between p zz and j p in a dilute 
vapor is the foundation of Eq. ©, in which the mass 
flux j m can be determined from pii q since pu q equals p zz 
in the vapor phase. As shown in Figs. [5l [H and [TOl the 
agreement between j m /m and the evaporation rate je is 
very good for all systems evaporating into a vacuum. 

For systems evaporating at small controlled rates, the 
vapor density is comparable to its equilibrium value and 
it is not expected that j p is equal to p zz even in the 
vapor phase, though p zz is still equal to py lq . As a result, 
the mass flux j m from Eq. @ is generally far from the 
evaporation rate mjE- A simulation of LJ monomers 
evaporating at Je = 1-0 x 10 _3 r _1 cr~ 2 and T& — 0.9e/£;B 
shows that j m /m from Eq. © is only about 50% of the 
measured je ■ 

V. COMPARISON TO KINETIC THEORY 

At the liquid/ vapor interface, the molecular exchange 
between the two phases occurs continuously in the form 
of evaporation and condensation. Results discussed 



above show that the molecular composition affects the 
evaporation rate to a large extent. It remains an interest- 
ing question if the change in molecular composition also 
affects the atomic kinetics during evaporation and con- 
densation, such as velocity distributions of the molecules. 
In equilibrium, all molecules leaving the interface to the 
gas phase or arriving at the interface from the gas phase 
have a normal velocity v z that satisfies the MB distribu- 
tion with a probability density function 



P(Vz) 



Mv z 



exp 



Mvl 
2k B T, 



(11) 



where M is the molecular mass and Ti is the temperature 
at the interface. Note that Tj equals Th in equilibrium but 
is generally lower than Th when there is net evaporation. 
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FIG. 12: (Color Online) Velocity distributions of v z for all 
molecules that leave (circles) and arrive (squares) at the in- 
terface for monomer systems at T;, = 0.9e/fcB under var- 
ious evaporation conditions: (a) liquid/vapor equilibrium; 
(b) controlled rate je = 2.5 x W~ 4 r~ 1 a~ 2 ; (c) controlled 
rate ]e = 1-0 x 10~ 3 r~ 1 <r _2 ; (d) in contact with a vacuum. 
Solid lines in (a) and (b) represent the MB distribution with 
T = T b = 0.9e/fc B . Two solid lines in (c) represent the MB 
distribution with T = 0.871e/fcB for molecules that leave the 
interface and Ti = 0.854e/fcB for molecules that arrive at the 
interface. Two solid lines in (d) represent the MB distribution 
with T — 0.744e/A:B for molecules that leave the interface and 
T = 0.292e/fcB for molecules that arrive at the interface. 

To analyze the velocity distributions, we consider 
a plane just outside the liquid/vapor transition layer. 
All molecules crossing this plane and entering (leaving) 
the transition layer are counted as incoming (outgoing) 
molecules. Their normal velocity v z is measured and the 
corresponding distribution function p(v z ) is calculated. 
The results for L J monomers are shown in Fig. [T^l Those 
for dimers and trimers are very similar and are not shown. 
The results are not sensitive to the position of the plane 
that molecules cross as long as it is outside but not far, 
within several <r, from the transition layer. 
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As expected, Fig. [T27a) shows that when the liq- 
uid/vapor interface is in equilibrium the MB distribution 
is satisfied by molecules that enter or exit the transition 
layer with a Ti that is the same as the bulk tempera- 
ture in the liquid and vapor phases. This is also the case 
when the evaporation rate is as low as 2.5 x 10 _4 r _1 (T -2 
[Fig. H2Tb)]. When the evaporation rate is increased to 
1.0 x lO-V-V" 2 [Fig. [12(c)], the MB distribution still 
holds, but corresponds to a T, that is slightly lower than 
Tf,, which is fixed by a thermostat near the lower con- 
fining wall. Furthermore, in this case slightly different 
Tj's have to be used for molecules that enter or exit the 
transition layer. 

The situation is quite different in the case of evap- 
oration into a vacuum, where the evaporation rate is 
very high. Figure IT2Td) shows that the velocity distri- 
butions of molecules entering or exiting the transition 
layer have very different TVs. For the outgoing molecules, 
Ti = 0.744e/fcB is consistent with the temperature of the 
liquid/vapor interface, which is lower than the bulk tem- 
perature Tb = 0.9e/fce of liquid far from the interface. 
This reduction of the temperature at the liquid/vapor 
interface is due to evaporative cooling. For the arriving 
molecules we found 7$ — 0.292e/fcB when p(v z ) was fit 
to a MB distribution. This temperature is much lower 
than that for the outgoing molecules and is only about 
1/3 of Tb. The strong asymmetry between Tj's for the 
outgoing and incoming molecules can be understood as 
follows. When the liquid evaporates into a vacuum, the 
vapor density is extremely low and it is very unlikely 
that a vapor molecule undergoes a collision with other 
molecules in the vapor phase and is reflected back to the 
liquid/ vapor interface. All molecules leaving the inter- 
face tend to move ballistically toward the deletion zone 
and get removed. This is particularly the case for out- 
going molecules with high velocities v z normal to the in- 
terface because they would require multiple collisions to 
be reflected, which is very rare. Therefore, most reflected 
molecules are those moving slowly. This leads to a strong 
imbalance in the distributions of v z of the incoming and 
outgoing molecules, which is manifest in the fact that 
the incoming molecules have a much lower Ti than that 
of the outgoing molecules. 

Results in Fig. H2T d) indicate that if the full velocity 
distribution is plotted for all the molecules approaching 
or leaving the liquid/vapor interface, one would expect 
some deviation from the MB distribution that is symmet- 
ric between v z and — v z . That is, the contribution from 
positive values of v z have more weight than that on the 
negative side. This deviation was also observed in pre- 
vious MD simulations of evaporation into a vacuum . 15 ! 33 
Our results further indicate that this asymmetry also oc- 
curs when the evaporation rate is high enough that the 
vapor density is substantially reduced from its equilib- 
rium value. 

Not all molecules arriving at the liquid/vapor inter- 
face from the gas phase condense into the liquid phase. 
A fraction is reflected back into the vapor phase. The 



condensation coefficient defined in Introduction quanti- 
fies this effect. It essentially represents the fraction of 
incoming molecules that indeed become attached to the 
liquid phase. It has been argued that the condensation 
coefficient is close to unity for monatomic liquids such as 
liquid metals^ and less than unity for polyatomic liq- 
uids because rotational motion of polyatomic molecules 
in the liquid state make it more difficult to accommodate 
newly arriving molecules^ 

The reflected molecules obviously contribute to the 
evaporation flux. Equivalently, not all evaporated 
molecules come directly from the liquid phase. The evap- 
oration coefficient is the ratio between the molecular 
flux due to real evaporation, i.e., from those molecules 
transformed into vapor from the liquid phase, and the 
maximum flux j max calculated from the MB distribu- 
tion. It is easy to derive the HK formula, j ma x = \nv, 
where n is the density of the saturated vapor correspond- 
ing to Ti and v is the mean molecular velocity. For 
liquid/vapor equilibrium the condensation and evapora- 
tion coefficients must be the same since the velocity of 
molecules arriving at the interface also satisfies the MB 
distribution. 

Hereafter, we consider the total evaporation flux as 
the sum of that due to the true evaporation of liquid 
molecules and another due to reflection of incoming vapor 
molecules. Results in Fig. [T27d) show that when a liquid 
evaporates into a vacuum, the reflection flux is greatly 
suppressed because of the depletion of vapor. In this 
case the total evaporation flux is very close to the true 
evaporation flux. This is the reason that Ishiyama et al. 
emphasized that the condensation coefficient can be de- 
termined without any ambiguity by measuring the spon- 
taneous evaporation flux from the simulations of evapo- 
ration into a vacuum^ 

Through MD simulations of LJ monomers, Tsu- 
ruta et al. discovered that the condensation coeffi- 
cient is actually an average quantityi 3 -^— For incom- 
ing molecules with different normal translational energy 
Ein,z = Mv z /2, their tendency to condensate into the 
liquid phase is generally different. A condensation prob- 
ability a c was defined to quantify this tendency and they 
suggested that a c depends on the normal velocity of in- 
cident molecules in the following functional form, 



1 — /3exp 



knT 



(12) 



where a and (3 are two parameters that depend on the 
properties of liquid and temperature. In this framework, 
the condensation coefficient discussed above is the av- 
erage of a c over the velocity distribution of incoming 
molecules and is denoted as a^. Since the velocity of all 
incoming molecules satisfies a MB distribution, aZ can 
be expressed as 



J +O ° a c (v z )p(v z )dv z 
a(l-0/2). 



(13) 
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FIG. 13: (Color Online) Condensation probability a c as a 
function of -Ein.z, which is the normal component of the trans- 
lational energy of incoming molecules. Data presented here 
are measured on liquid/ vapor equilibrium systems at various 
temperatures. However, results are similar if the evapora- 
tion rate is fixed at a small value and the liquid/vapor inter- 
face is near equilibrium, (a) monomers: Tj, = OSe/ku (tri- 
angles), 0.9e/A)B (squares), and l.Oe/ks (circles); (b) dimers: 
Tb = 1.2e/&B (triangles), 1.3e/&B (squares), and 1.4e/fcB (cir- 
cles); (c) trimers: Tt = lAe/ks (triangles), 1.5e/fcB (squares), 
and 1.6e/feB (circles). In this plot a conversion factor e/fcs = 
119.8K (for argon) is used. 



In our simulations we grouped incoming molecules by 
the normal component v z of their velocity and measured 
their condensation probability a c . However, an ambi- 
guity arises about what incoming molecules should be 
counted as condensed ones since all molecules arriving 
at the transition layer will eventually leave after some 
period of time. We counted those incoming molecules 
that stay in the liquid phase or the liquid/vapor transi- 
tion zone for at least At = 25r. Since a typical value 
r is approximately 2ps, 56 At is about 50ps. This choice 
of At is somewhat ad hoc. However, previous MD sim- 
ulations have established that the characteristic time of 
energy excitation during evaporation or relaxation dur- 
ing condensation is roug hly 50 - 70ps.2£ For this reason 
our choice of condensation time scale as 50ps should be 
reasonable. This choice is further substantiated by the 
fact that the measured a c does not change significantly 



even if At is increased by a factor of 2 to 50t. All data 
presented here are for At = 25r. After a c was obtained, 
the a c ~ v z was fit to Eq. (TT2"|) with a and j3 as fitting 
parameters, from which the condensation coefficient 
was calculated using Eq. (fl~3|). 

Results for a c (v z ) are shown in Fig. [T3] It indicates 
that a c indeed follows the functional form in Eq. (|12l) . 
regardless of the change in molecular composition from 
monomers to dimers and to trimers. Vapor molecules 
arriving at the liquid/ vapor interface with a small v z are 
less likely to condense and join the liquid state. This can 
be understood by a simple physical picture. When the 
arriving molecules collide with liquid molecules, they are 
more likely to be reflected after one or several collisions 
if their velocity is small. Arriving molecules with a large 
velocity v z are able to penetrate the liquid region deeper, 
and have a bigger chance to survive multiple collisions 
and to remain in the liquid phase for at least a period 
-At. 

The condensation probability a c being a function of v z 
has another implication. Since all molecules leaving the 
liquid/ vapor interface are either those that truly evap- 
orate from the liquid or those that have been reflected, 
and the velocity of all leaving molecules satisfies the MB 
distribution (as shown in Fig. H"2"1) . the normal velocity 
distributions of the truly evaporated molecules and re- 
flected ones are modified from the MB distribution. Par- 
ticularly, the MB distribution is modified by the conden- 
sation probability a c for the truly evaporated molecules 
and by 1 — a c for the reflected ones. Thus they have the 
following forms after a proper normalization, 30 
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(15) 

where the subscripts e and r stand for "truly evaporated" 
and "reflected", respectively. 

The distribution functions p e (v z ) and p r (v z ) for the 
truly evaporated and reflected molecules were measured 
directly in simulations and provided another way to de- 
termine a and /3, and eventually the condensation coeffi- 
cient ~a~ c . One example is shown in Fig.[l4j which includes 
data on p(v z ), p e (v z ), and p r (v z ) of monomers evaporat- 
ing at a small fixed rate. The successful fits to the MB 
and modified MB distributions in Eqs. (jlll) . (fl4)) . and 
(|15|) indicate that the condensation probability a c pro- 
vides a reasonable quantitative measure of the conden- 
sation and evaporation processes. Furthermore, we com- 
pared the condensation coefficient aZ determined either 
using a c (v z ) or using p e {v z ) and p r (v z ). For systems in 
equilibrium or not far from equilibrium such as evapora- 
tion with small controlled rates, the two results generally 
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FIG. 14: (Color Online) Velocity distributions of v z of all 
molecules which left (triangles), were evaporated (squares), 
or were reflected (circles) from the interface for LJ monomers 
at Tb = 0.9e/ltB evaporating at a rate ]e = 2.5 x 10 _4 r -1 <7~ 2 . 
Solid lines are fits to the MB distribution p(v z ) in Eq. (|ll|l and 
the modified MB distributions Pe(v z ) and p r (v z ) in Eqs. (|14p 
and (|15p . respectively. 



agree, which further validates the kinetic model of evapo- 
ration described above. For systems far from equilibrium 
such as evaporation into a vacuum, it is generally very 
difficult to measure aZ in these ways because there are 
essentially no incoming molecules and thus no reflected 
flux. In this case, it might be possible to determine aZ 
by directly measuring the evaporation flux and compar- 
ing to the theoretical maximum flux j max , as suggested 
in Ref. [H. However, the temperature of the liquid/vapor 
interface needs to be determined first in order to calcu- 
late j max . This is not an easy task as the interface is 
moving during evaporation and the temperature is sensi- 
tive to the location of measurement in the region around 
the interface. 

Figure IT51 shows results on the condensation coefficient 
ctZ for LJ monomers, dimers, and trimers. The agreement 
between results from the condensation probability and 
the velocity distributions is clear. Generally, aZ decreases 
with a increasing T. For monomers in temperature range 
we could study, the reduction is approximately linear. 
The expected crossover in aZ to unity for lower tempera- 
tures is not observed for monomers since for lower T the 
system crystallizes. But the crossover can be identified 
in the data from Ref. |H, which used the Dymond- Alder 
potential for argon in the low temperature range. For 
dimers and trimers, the expected crossover of ~ T at 
lower temperatures is evident. 

Figure [TS] shows that the value depends on the evap- 
oration condition. This is obvious if we examine the 
monomer data at T h = 0.9e/fce. The value of aZ in liq- 
uid/vapor equilibrium is slightly lower than its value at 
the evaporation rate je = 2.5 x lCPV^er -2 , and both 
are clearly lower that the value at the larger evaporation 
rate ]e = 1.0 x 1Q _3 t _1 (7 _2 . However, if we take into 
account the fact the temperature Tj of the liquid/ vapor 
interface get lower when evaporation is stronger and use 




0.4 
0.2 
0.0 



3.0 0.1 0.2 0.3 0.4 0.5 0.6 



4 0.5 0.6 0.7 0.8 0.9 

T/T c 



1.0 



FIG. 15: (Color Online) Condensation coefficient vs. tem- 
perature for monomers, dimers, and trimers under various 
evaporation conditions. For monomers: circles (O) an d tri- 
angles (A and v) are measured using a c (v z ); open squares 
(□), pluses (+), and crosses x are measured using p e (v z ) 
and Pr(v z )\ circles (O) an d open squares (□) are for liq- 
uid/vapor equilibrium systems at T = T b — 0.8e/fcB, 0.9e/fcB, 
and 1.0e/fcB; upward triangles (A) and pluses (+) are for 
evaporation rate je = 2.5 x 10~ 4 T~ 1 a~ 2 at T, = 0.9e/£;B 
and 1.0e/fcB; downward triangles (v) and crosses (x) are for 
evaporation rate js = 1.0 x 10 -3 t - V" 2 at T& = 0.9e/fcB- For 
data measured during evaporation, the actual temperature T 
of the liquid/vapor interface is used. The data for dimers 
and trimers are shown with solid squares (■) and solid tri- 
angles (A), respectively, and both are calculated using a c (v z ) 
from liquid/vapor equilibrium systems. The dashed line is 
a linear fit to the monomer data at high temperature. The 
dotted lines are guides to the eye. Our data of LJ monomers 
agree with those from Ref. [3g [shown with open diamonds 
(0)], which were determined by directly measuring the evap- 
oration flux into a vacuum. The inset plots aZ as a function of 
the translational length ratio, (Vi/K,) 1 ^ 3 , between the vapor 
and liquid phase. The solid line in the inset is the prediction 
of the transition state theory (see Ref. l3lT) . 



Tj instead of T fc for the horizontal axis, all data fall nicely 
onto the same master curve indicated by the straight 
line in Fig. 1151 The reason that is larger under 
stronger evaporation is that the interface temperature 
decreases due to evaporative cooling and lower temper- 
atures corresponds to larger a^. This trend is also in- 
dicated by the data at T h = l.Oe/fce for an equilibrium 
interface (thus Tj = T h ) and for a fixed evaporation rate 
j E = 2.5 x lO-V-V- 2 (thus T, < Tj,). In the latter W c 
is found larger. 

Our results show that the condensation coefficient aZ 
is substantially less than unity for LJ monomers above 
the triple-point temperature, and is higher for LJ dimers 
and trimers at the same reduced temperature T/T c . This 
is in contrast to the traditional viewpoint that the con- 
densation coefficients decreases when the molecular com- 
position changes from monatomic to polyatomic. In this 
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viewpoint, the rotational motion of polyatomic molecules 
acts as a constraint to make the vapor condensation 
less likely. However, the short chain molecules we have 
simulated are relatively flexible. When two polyatomic 
molecules come in contact with the interface, their at- 
tractive interaction is stronger than that between two 
monatomic molecules because of the increase in the num- 
ber of interacting sites. Furthermore, because of the flex- 
ibility, rotational degrees of freedom play a less impor- 
tant role. This leads to stronger cohesion for polyatomic 
chain molecules. It is therefore not surprising that 
is higher for LJ dimers and trimers at the same T/T c 
because it basically represents how likely a molecule can 
bind to its liquid phase. Tsuruta and Nagayama also 
found that translational degrees of freedom dominate in 
the evaporation and condensation of water^i 

Using transition state theory, Nagayama and Tsuruta 
derived a theoretical expression for in the case that the 
evaporation and condensation processes are dominated 
by translational motion of molecules^ In this frame- 
work, has a universal functional dependence on a 
translational length ratio, which is defined by the cubic 
root of the free volume ratio between the liquid and vapor 
phase, (Vi/Vg) 1 / 3 . The results for W c versus (Vz/Vg) 1 / 3 
are shown in the inset of Fig. [15] All data collapse onto a 
master curve except for dimers at low temperature. How- 
ever all data are above the theoretical line. Most data 
for in Ref. [3l| are also above the predicted curve for 

or c vs. (yi/Vg) 1 ' 3 . 

VI. SUMMARY AND CONCLUSIONS 

We have used MD simulations to investigate fluids 
made of LJ monomers, dimers, and trimers. The phase 
diagrams for the dimer and trimer systems were deter- 
mined and their evaporation processes were simulated. 
Our results show that a simple change from monomer to 
dimer or trimer molecules has a strong effect on the evap- 
oration rates. The evaporation rate of monomers when in 
contact with a vacuum is extremely high and there are 
strong evaporative cooling and liquid density enhance- 
ment near the liquid/vapor interface. All these effects 
are greatly reduced in dimer and trimer fluids. A phys- 
ical explanation is provided on the basis that cohesion 
gets stronger for fluids made of longer chain molecules 
because of the increase in the number of interacting sites 
between two close molecules. The flexibility of the linear 
chain also suppresses the importance of rotational motion 
of molecules, which tends to reduce the cohesion. The 
net outcome of the competition between these factors is 
that the surface tension is higher and the evaporation 
slows down at the liquid/vapor interface for fluids made 
of longer chains. 

The measured evaporation rates for liquids evaporat- 
ing into a vacuum were compared to the modified HK for- 
mula derived by Holyst and Litniewski. Good agreement 
was found for not only monatomic liquids composed of 
monomers, but also molecular liquids composed of dimers 



and trimers. It was confirmed that mechanical equilib- 
rium, i.e., a constant component of the stagnation pres- 
sure normal to the liquid/vapor interface, holds for all 
systems at all times during evaporation. Furthermore, 
for liquids evaporating into vacuum, the momentum flux 
contribution dominates in the stagnation pressure in the 
vapor phase, which was first shown by Holyst and Lit- 
niewski to lead to the modified HK formula^ 

Because of evaporative cooling, the temperature of the 
liquid/vapor interface is lower than the bulk temperature 
of liquid. It is further observed that the temperature in- 
creases in the vapor phase with the distance from the liq- 
uid/vapor interface. This makes the interface the coolest 
point in the system. The phenomenon that the vapor 
phase close to the interface has a temperature higher 
than that at the interface is qualitatively consistent with 
a previous experimental measurement on water 

We measured velocity distributions of molecules at 
the liquid/vapor interface. As expected, in liquid/vapor 
equilibrium they follow the MB distribution. When the 
evaporation rate is controlled at a small value, the distri- 
bution also has a MB form but may correspond to a tem- 
perature that is slightly lower than the bulk temperature 
of either liquid or vapor phase. When a liquid evaporates 
into a vacuum, the molecules arriving at the interface 
have a velocity distribution that corresponds to a tem- 
perature much lower than that for the molecules leaving 
the interface. Furthermore, both temperatures are lower 
than the bulk temperature of the liquid phase far from 
the interface. In this case, evaporated molecules move 
away balistically in the vapor phase without encounter- 
ing any impedance because the vapor density is extremely 
low. It is not surprising that local thermal equilibrium 
is not reached in this extreme situation. However, our 
data also showed that under moderate evaporation rates, 
which are usually the case in experiments, the vapor den- 
sity remains comparable to its equilibrium value and the 
hypothesis of local thermal equilibrium becomes valid. 

The condensation coefficient 3£ was determined by 
measuring the probability of reflection of molecules ar- 
riving at the liquid/ vapor interface from the vapor phase. 
This probability generally depends on the normal com- 
ponent of the velocity of the arriving molecules. The 
functional form was confirmed to be exponential as in 
Eq. (|1 lj) . consistent with the model of condensation prob- 
ability of Tsuruta et al.M-~— Fitting the probability to 
Eq. (|lip gave two parameters a and j3, from which ~ai 
was calculated. 

The apparent evaporation flux contains two contribu- 
tions. One is from molecules that truly evaporate from 
the liquid phase, and another is from molecules arriving 
at the interface from the vapor phase but being reflected 
back. The velocity distributions of these two groups were 
measured and compared to modified MB distributions 
that depend on a and /3. From these velocity distribu- 
tions, a, /? and eventually the condensation coefficient 
were also determined. 

The condensation coefficients aZ measured with the 
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above two methods are very close, which serves as a vali- 
dation of the kinetic model of Tsuruta et al.~~— We also 
found that decreases with increasing temperature T, 
consistent with the intuition that condensation becomes 
less likely at higher T. At the same reduced temperature 
T/T c , trimer fluids have the largest and monomers 
have the smallest. This is in contrast to the traditional 
viewpoint that monatomic liquids have a condensation 
coefficient close to unity and polyatomic liquids have a 
value less than unity. However since the chain molecules 
simulated here are very flexible, we do not think our 
results are in contradiction this viewpoint. The chain 
flexibility makes less prominent the rotational degrees of 
freedom that tend to frustrate the accommodation of in- 
coming molecules to the liquid phase. Increasing in chain 
length leads to a stronger cohesive interaction between 
chain molecules. As a result, the condensation coeffi- 
cient is higher for liquids made of molecules of longer 
chains. However, in most low molecular weight liquids, 
the molecules are relatively rigid, and the rotational mo- 
tion plays a more important role in determining liquid 
cohesion and molecular orientation at the liquid/ vapor 
interface. This can reduce the condensation coefficient, 
as suggested by some experiments. It remains an inter- 
esting open question to see if our observation that larger 
~aZ occurs for longer chains can be borne out by study- 



ing liquids made of a series of chain molecules such as 
hydrocarbons. 

The condensation coefficients for monomers, 

dimers, and trimers almost collapsed onto a single master 
curve when plotted against a translational length ratio, 
(Vi/Vg) 1 ^ 3 , between the vapor and liquid phase. How- 
ever, all of our data are consistently above the theoretical 
prediction of Nagayama and Tsuruta based on transition 
state theory^ Note that Vi/V g is just the inverse of the 
density ratio. This deviation indicates that the struc- 
ture of a liquid/ vapor interface is more complicated than 
that of a simple material-dividing plane characterized by 
a single parameter, (Vi/Vg) 1 ^ 3 . 
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